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Abstract 

We study a simple lattice model of shear-induced clustering in two dimensions in 
which clusters of particles aggregate under an imposed shear flow and fragment 
stochastically. Two non-equilibrium steady states are identified: an unjammed state 
and a jammed state characterised by a system-spanning cluster. A discontinuous 
jamming transition with strong hysteresis occurs as the shear rate is increased or 
fragmentation rate decreased. We study the kinetics of jamming and measure power 
law cluster size distributions. We also consider some general simulation issues in- 
cluding the role of Galilean invariance. 



1 Introduction 



When subjected to shear flow, many complex fluids exhibit shear thickening 
behaviour where viscosity increases with shear rate. Shear thickening is often 
discontinuous and may be accompanied by hysteresis and structural changes. 
Systems where such behaviour has been observed include concentrated col- 
loidal suspensions [1] and rodlike micelle solutions [2]. 

In the case of colloidal suspensions, it has been suggested that shear thickening 
is the result of the formation of hydrodynamic clusters [3]. These clusters 
comprise large numbers of particles bound together by short-range lubrication 
forces and form along the compression axis of the shear. Indeed a detailed 
analysis of the dynamics of the aggregation of such rotating clusters, in the 
absence of Brownian motion, has recently been performed [4] and suggests 
that the discontinuous jump in viscosity is due to a jamming transition where, 
essentially, a log-jam of clusters occurs. 

Jamming transitions have also been studied in an apparently different context, 
that of driven diffusive systems [5]. In one dimensional models, jams have 
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been induced by the presence of defects or disorder [6] or by a non-conserved 
quantity mediating an effective long-range interaction [7]. In two-dimensional 
(2d) driven lattice gases, gridlock effects have been studied [8,9]. 

Recently, Santra and Herrmann introduced a simple 2d lattice model of gela- 
tion under shear [10], which provides a link between colloidal systems and 
driven lattice gases. They showed that within the model one finds dynamical 
scaling as clusters of particles aggregate irreversibly, and that the clusters are 
aligned roughly along the compression axis of the shear. The model has been 
extended to take into account the effect of rotation of the clusters [11]. 

In the present work, we extend the model of [10] in a different direction by 
allowing clusters to fragment stochastically. This models thermal fluctuations 
in a simple manner and gives rise to the possibility of non-trivial steady states, 
in contrast to the absorbing state consisting of a single cluster which is always 
attained with the models of [10] and [11]. Indeed, as we shall show below, we 
find two types of steady state - unjammed, where clusters have some typical 
size much smaller than the transverse size of the system, and jammed, in which 
there are clusters which span the system. We provide evidence of a discontin- 
uous transition between the two states, with hysteresis, as the fragmentation 
rate or the shear rate is varied. We study the kinetics of the jamming process 
and find evidence for dynamical scaling for both zero and non-zero fragmenta- 
tion rate. Although there is an intriguing similarity with the colloidal jamming 
problem studied in [4], in practice the transitions we observe arise at very low 
concentration. In this respect the behaviour is more like that of the micellar 
problem of [2] . 

The paper is organised as follows. In section 2 we define the model that we 
study and simulation results are presented in section 3. In section 4 more 
general issues concerning the simulations are raised. We conclude in section 5. 



2 Description of the Model 



The model is defined on a 2d lattice of length L and width W, illustrated 
by an example configuration in figure 1. Each lattice site is either occupied 
by a particle or empty. The boundary perpendicular to the direction of flow 
is periodic but no particles may pass through the boundaries parallel to the 
direction of flow - they are "hard walls" . A simple shear flow is imposed by 
assigning a velocity v (y) to each "lane" in the direction of flow given by 

v(y) = v(0) + iy (1) 
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Fig. 1. Example lattice configuration. The shear flow is represented by the arrows. 

where y runs from to W — 1. Each particle in lane y is said to have velocity 
v(y). Hence, the shear rate is 
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W - 1 



(2) 



where i> m i n = v(0) and t> max = t>(W — 1). We consider only non-negative veloc- 
ities from now on. There are N particles in the system and the concentration 
is 



We define a cluster to be a set of particles connected successively by nearest- 
neighbour bonds on the lattice. Therefore, the nearest neighbour sites belong- 
ing to an individual particle are either empty or occupied by particles in the 
same cluster. The velocity u of a cluster is defined as the mean of the velocities 
of its constituent particles. The proximity of a wall does not affect the velocity 
of a cluster so that one may think of the particles "slipping" relative to the 
walls. (In contrast, in order to set up a shear flow in a fluid by moving one 
or both of the walls, a no-slip condition between the fluid and the walls is 
required.) 

A (continuous) "counter" variable q is associated with each cluster. In each 
simulation time-step, three steps (flow, cluster merging and fragmentation) 
are performed in the following sequence: 

Flow The counter q for each cluster is incremented by u, the cluster velocity. 
If q > 1 for any cluster, all the particles in that cluster move one site forward 
(i.e. in the direction of flow) and q — ► q — 1. If q < 1, no action is taken. 

Cluster merging At the end of the previous step, some clusters may have 
become connected and these clusters are now merged. In merging two clus- 
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ters, the resulting single cluster is assigned a counter which is the weighted 
average of the counters of the two constituent clusters labelled i and j: 



_ ni q t + n jqj 

where rii is the number of particles in cluster i. We have checked that this 
ensures that the average velocity of clusters in lane y is indeed v(y). 
Fragmentation This random sequential step allows single particles to leave 
clusters, possibly resulting in the break up of clusters. The following update 
rules are repeated N times: 

(1) A particle is chosen at random. 

(2) A direction for the particle to move is chosen at random from the four 
possibilities. 

(3) If the particle is not blocked by another particle in the chosen direction, 
it moves one site in that direction with probability T, where T may be 
thought of as the "temperature" of the system. 

(4) The moved particle becomes a single-particle cluster with counter vari- 
able equal to that of the cluster it previously belonged to. If it adjoins 
any other clusters, this single-particle cluster is then merged as detailed 
previously. 

The fragmentation update of one particle may result in one of five outcomes: 
no change, a particle leaving its cluster, a particle moving to a different 
location in the same cluster, movement of a single-particle cluster or a single- 
particle cluster joining another cluster. 

The three physical parameters present in the model are concentration c, shear 
rate 7 and "temperature" T. 

Simple shear flow is Galilean invariant, by which we mean that adding a 
constant velocity at every point in the system has no effect on the physics. It 
is therefore desirable that any model of simple shear is as close as possible to 
being Galilean invariant. We discuss this issue in section 4 where we show that 
lattice effects break Galilean invariance in the model. In addition, we remark 
that the model is translationally invariant only in the direction of flow. The 
fixed walls break translational invariance normal to the direction of flow. 



3 Results 



3.1 Steady State Behaviour 



From computer simulation of the model defined in section 2, we have found 
two types of steady state. The first is an unjammed state, where many small 
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Fig. 2. Snapshot of an unjammed 1200x200 system in the steady state. The simu- 
lation parameters were c = 0.07, T = 0.0045 and 7 = 5/1990 ~ 0.0025. The flow is 
from left to right and the velocity increases from bottom to top. 



Fig. 3. Snapshot of a jammed 1200x200 system in the steady state. The sim- 
ulation parameters were the same as for figure 2 with the exception that here 
7 = 7/1990 ~ 0.0035. 

clusters are evenly distributed throughout the system. A snapshot of a typical 
unjammed system is shown in figure 2. There are some clusters which are 
elongated perpendicular to the direction of flow but the typical cluster size is 
very small. 

A second type of steady state exhibited by the model is characterised by the 
presence of clusters which span the system; we refer to it as "jammed" state 
(although the spanning cluster is not stationary). These spanning clusters 
contain a finite fraction of the particles in the system. A snapshot of a typical 
jammed system is shown in figure 3; the system is strongly inhomogeneous in a 
number of respects. Since the spanning cluster is moving with a velocity close 
to that of the particles in the centre of the channel, the rate of aggregation is 
least in the centre and greatest at the walls - near the upper wall, small clusters 
move more quickly than the spanning cluster and join it from behind whereas 
near the lower wall, small clusters move more slowly than the spanning cluster 
and aggregate from the front. Hence the spanning cluster is aligned at an angle 
to the direction of flow and decreases in thickness from the edges to the centre 
of the channel. Accordingly, the density in the rest of the system is greatest in 
the centre and least at the edges. Some quite large clusters are visible close to 
the spanning cluster - these are probably pieces of the spanning cluster which 
broke off some short time previously. The somewhat unusual, inhomogeneous 
nature of the jammed steady state in the finite systems we study is primarily 
a result of having hard wall boundaries parallel to the direction of flow rather 
than any effect of the underlying square lattice. 

In order to discuss the behaviour of the model more quantitatively, we intro- 
duce the cluster size distribution n(m,t), the number of clusters containing 
m particles at time t. The moments M k (t) of the cluster size distribution are 
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Fig. 4. Plot of s against 7 for a 900x150 system with c = 0.07 and T = 0.0045. 
The direction of the hysteresis loop is shown by the arrows. 

given by 

00 

M k (t)= J2m k n(m,t) (5) 

m=l 



so that Mi = N, the number of particles in the system and M is the number 
of clusters in the system. We define 

M k+1 (t) 



as measures of typical cluster size [12] for k > 0. In discussing steady state 
properties of the system, we use the average cluster size s as our measure 
of typical cluster size since this quantity is not dominated by the spanning 
clusters in jammed systems. In contrast, for k > 0, s k is of the order of the 
size of a spanning cluster in the jammed state. We use si when discussing the 
kinetics of jamming in section 3.2. 

We now provide evidence of a strong discontinuous phase transition between 
jammed and unjammed states as the parameters T and 7 are varied. Figure 
4 shows a plot of s against 7. Strong hysteresis is apparent, with the direc- 
tion of the hysteresis loop shown by the arrows. As 7 is increased from zero, 
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Fig. 5. Log-log plot of so against T for a 900x150 system with c = 0.06 and 
7 = 6/1490 ~ 0.004. The direction of the hysteresis loop is shown by the arrows. 

So increases smoothly at first but jumps discontinuously to a larger value, 
corresponding to a jammed system. When 7 is then decreased, the jammed 
state remains stable (for at least 2 x 10 5 time-steps) for a range of 7 val- 
ues in which an unjammed state is also stable (for at least 10 5 time-steps). 
The system becomes unjammed again for sufficiently small 7 {i.e. an initially 
jammed configuration becomes unstable). One can also see from figure 4 that, 
for jammed systems, so increases roughly linearly with increasing shear rate. 

Similar behaviour is found when T is varied with 7 and c held fixed. Figure 
5 shows a log-log plot of sq against T. As in figure 4, strong hysteresis is 
apparent. There is a large range of values of T for which both unjammed and 
jammed states are stable (on a time-scale of at least 2 x 10 5 time-steps for 
jammed systems and at least 10 5 time-steps for unjammed systems.) 

Figure 6 shows a log-log plot of the steady state cluster size distribution n(m) 
for both unjammed and jammed systems at the same value of T = 0.006. In 
the homogenous system, n{m) decays very rapidly and the large m behaviour 
appears to be exponential (it decays faster than any power law). In the jammed 
system, the size distribution of small (non-spanning) clusters only is shown; 
the distribution as a whole is bimodal, with the second peak (not shown) 
corresponding to the spanning cluster. The figure suggests that for the small 
clusters, n{m) decays as a power of m (for m greater than about 10) with an 
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Fig. 6. Log-log plot of the steady state cluster size distribution for both unjammed 
and jammed systems with T = 0.006. The system size is 900x150, c = 0.06 and 
7 = 6/1490 ~ 0.004. The dashed line is proportional to m~ 2 - 75 and is shown to 
guide the eye. 

exponent close to —2.75. However, n{m) decays faster than a power law for 
m greater than about 100. It is not clear whether this is a finite system-size 
effect or the true asymptotic behaviour. For jammed systems, the situation is 
further complicated by the fact that the cluster size distribution is spatially 
inhomogeneous (see figure 3). 

We have found two different regimes of behaviour in the jammed state as T is 
varied. Figure 7 shows the variation of Sq with time in the (jammed) steady 
state with T = 0.007 (the other parameters are the same as for figure 5). 
The average cluster size s appears to fluctuate about a value close to 3.25 
for some time before increasing suddenly to more than 4.0; it then gradually 
falls back to about 3.25. We have found that this behaviour is due to splitting 
and reforming of the spanning cluster. We have already seen in figure 3 that a 
spanning cluster is narrowest at its centre. If T is sufficiently small, the span- 
ning cluster does not split. However, when T is sufficiently large the spanning 
cluster may split in two, whereupon the two resulting clusters separate. At 
this point, the two large clusters grow by aggregating with smaller clusters 
before meeting again and reforming the spanning cluster (due to the periodic 
boundary in the flow direction). Upon reforming, the spanning cluster is larger 
(it contains more particles) than when it split; this is because the two large 



8 




Fig. 7. Average cluster size sq against time for a 900x150 system with T = 0.007, 
c = 0.06 and 7 = 6/1490 ~ 0.004. 

clusters formed as a result of the split grow significantly before recombining. 
This explains the sudden rises in average cluster size seen in figure 7. The 
newly reformed spanning cluster gradually loses particles before splitting once 
more. The steady state is cyclic - a spanning cluster is stable only in the 
sense that once it splits, it always reforms. Close examination of figure 5 re- 
veals a kink in the jammed branch of Sq(T) around T = 0.007. This is caused 
by the change in behaviour from a non-cyclic steady state to a cyclic one. 
If T is further increased beyond about 0.018, the two large clusters formed 
when a spanning cluster splits may become unstable and break up (due to 
fragmentation dominating shear-induced aggregation). An unjammed system 
then results and there is no longer a jammed steady state. It is not clear to 
what extent the splitting of spanning clusters is affected by system size. 



3.2 Jamming Kinetics 



The above concludes our discussion of the steady state properties of the model; 
we now consider the approach to the steady state when the system sponta- 
neously jams from an initial condition of particles positioned at random. Fig- 
ure 8 shows a snapshot of a system in the process of jamming 7000 time-steps 
after the particles were positioned at random. Clusters with a large range of 
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Fig. 8. Snapshot of a 1200x200 system undergoing jamming at t = 7000. The 
particles were positioned at random at t = 0. The simulation parameters were 
c = 0.06, T = 0.002 and 7 = 8/1990 ~ 0.004. 
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Fig. 9. Plot of igei against T for various system sizes. The simulations had c = 0.06 
and 7 = 8/1990 ~ 0.004. The particles were positioned at random at t = 0. We 
averaged over betwen 10 and 20 independent runs. 

sizes are visible. 

We define the gel time t ge \ as the time at which a spanning cluster first appears 
in the system. Santra and Herrmann [10] found (by extrapolation) that, for 
T = 0, tg e i is finite for an infinite system. Figure 9 shows a plot of t ge \ against 
T for different system sizes. It appears that as system size is increased for 
the fixed aspect ratio L/W = 10, t ge \(T) converges. A fit to the data for 
the 2000 x 200 system indicated that t gc \ increases more quickly than exp(T) 
and it is possible that it diverges for some finite value of T; this is consistent 
with the presence of a phase transition at finite T. However, due to the large 
amount of computational effort required to obtain t ge \(T), we have been unable 
to confirm this unambiguously. Also, figure 9 shows that t ge \ depends on the 
aspect ratio of the system since t ge \ for the 10000 x 200 system is consistently 
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Fig. 10. Log-log plots of n{m) at different times for T = 0.001. The system size 
was 10000x200 and c = 0.06, 7 = 8/1990 ~ 0.004. The particles were positioned at 
random at t = 0. We averaged over 12 independent runs. The data has been shifted 
vertically for clarity. 

smaller than for the 2000 x 200 system. 

In a gelling system, we expect the pre-gelation cluster size distribution to have 
the scaling form [12,13] 



for large cluster sizes, with the typical cluster size si(t), given by (6), diverging 
as t — > t ge \. Since the number of particles N is conserved, we must have r > 2. 
The above scaling ansatz is expected to apply for an infinite system. In our 
finite systems, we find evidence of such scaling for T > although finite-size 
effects are large. 

Figure 10 shows a log-log plot of n(m) for T = 0.001 at four different times 
(note that the curves have been shifted vertically for clarity). Power law be- 
haviour appears gradually and is present over the largest range of cluster sizes 
at t* ~ 4500. We define t* to be the time at which we see power law behaviour 
of n(m) over the maximum range of m. By t = 5000, spanning clusters have 
formed and a "bump" has begun to form in the tail of the distribution. By 



11 



fitting a power law to the n(m) for large m at t — t*, we find the exponent 
r ~ 1.9. Since for an infinite system r must be greater than 2, we conclude 
that finite-size effects are strong. Indeed, this is not surprising since we find, 
for T < 0.0025, that t* ~ t ge \ (as expected in gelling systems) so that there 
are clusters with size of order the system size present at time t*. However, for 
T = 0.003 and T = 0.0035 we find that t* is somewhat smaller than t ge \, an 
anomaly which we are unable to explain but which may be due to finite-size 
effects. 




Fig. 11. Log-log plots of n(m) for various values of T after a time t*(T). The system 
size was 10000x200 and c = 0.06, 7 = 8/1990 ~ 0.004. The values of t* are 3400, 
4500, 7000, 8500, 10000 and 12000 (in order of increasing T). We averaged over 
12 runs for each value of T. The data has been shifted vertically for clarity. The 
dashed line is proportional to mT 2 and the solid line is proportional to m~ 2 ' 8 ; they 
are shown to guide the eye. 

Figure 11 shows log-log plots of n(m) at t* for several values of T. A power 
law cluster size distribution (for m larger than about 10) is present for all 
values of T. The slopes of the power law regions decrease with increasing T 
(for T > 0), suggesting that the exponent r varies significantly with T and 
hence that it may be non-universal. Futher careful study of finite-size effects 
would be required to clarify this point. 
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4 Simulation Issues 



The "counter" scheme [10] for performing cluster flow described in section 2 
is not the only one that can be imagined. However, we have found that it has 
important advantages over other schemes. 

For example, a simpler way of implementing cluster motion in the direction 
of flow is, at each time-step, to let each cluster move forward one site with 
probability v, where v is the velocity of the cluster. However, there are seem- 
ingly insurmountable difficulties with this stochastic implementation of cluster 
flow. To illustrate this point, let us consider how clusters aggregate and grow. 
Aggregation occurs when one cluster catches up with a slower neighbouring 
cluster. We wish the imposed shear to be the dominant factor in the ag- 
gregation. However, if clusters move forward stochastically, as in the scheme 
just described, fluctuations (in addition to the imposed velocity gradient) will 
cause clusters to aggregate. For example, a cluster with velocity v\ may may 
move forward in several successive time-steps at the same time as a cluster 
with velocity v 2 > v i does not move at all, possibly resulting in the spurious 
aggregation of the clusters. Indeed, we have found that when the stochastic 
implementation of cluster flow just described is used (instead of the counter 
scheme described in section 2), stochastic aggregation of clusters dominates 
aggregation due to the imposed shear for practical parameter values. 

A further advantage of the counter dynamics for flow is that it is much closer to 
being Galilean invariant than a stochastic implementation of cluster flow. This 
is because, for stochastic flow, the size of the fluctuations in the motion of a 
cluster, and hence the rate of cluster-cluster aggregation, is strongly dependent 
on the absolute magnitude of the velocity v(y). It is well-known that the 
variance of the displacement of an independent random walker hopping in 
one direction with probability per time-step v is proportional to v(l — v); we 
therefore expect the rate of cluster-cluster aggregation for stochastic cluster 
flow to be greatest for v = 0.5 and smallest for v — and v — 1. We have 
confirmed that this is indeed the case in simulations. This is a potentially 
serious problem in simulation of a shear flow since v varies across the width 
of the lattice. 

Since the counter implementation of cluster flow is deterministic, the problems 
peculiar to stochastic cluster flow just discussed are circumvented. The counter 
scheme is, however, not without its limitations. Figure 12 shows the average 
cluster size of clusters whose centre of mass lies in lane y as a function of 
v(y) for two unjammed systems with different concentrations. The average 
cluster size is smallest for v — and v — 1 due to the presence of the fixed 
walls. However, the cluster size grows quite slowly as one moves in from either 
wall and we have checked that the reason for this is that v is in the vicinity 
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Fig. 12. Mean size of clusters having velocity v in two unjammed systems (averaged 
over 10 5 time-steps in the steady state). In both cases, the system size is 400 x 100 
with Vrai n = and f max = 1- The top curve is for c = 0.15 and T = 0.045 while for 
the bottom curve, c = 0.07 and T = 0.02. 

of either or 1. Similarly, a deep minimum in the average cluster size is 
evident around v = 1/2, with less pronounced minima also evident around v = 
1/3,2/3, 1/4,3/4, etc. Clearly, the counter updating scheme has a systematic 
effect on the cluster-cluster aggregation rate, breaking Galilean invariance by 
introducing a velocity dependence of the cluster-cluster aggregation rate. Not 
surprisingly, the effect diminishes as c is decreased. In our simulations, we 
have attempted to minimise the effect of the spatial inhomogeneity caused by 
the updating scheme by working at low concentration and avoiding velocities 
less than 0.1 or greater than 0.9. For the concentration values we have used 
(c < 0.07), it is clear that the shear is the dominant factor driving cluster- 
cluster aggregation since the average cluster size varies significantly with shear 
rate (see figure 4). 

It is clear that lattice models are not ideally suited to the simulation of a 
Galilean invariant phenomenon such as shear-induced aggregation, although 
they are a relatively simple and efficient approach to the problem. However, 
there are physical situations where Galilean invariance is not respected, an 
example being (multi-lane) traffic flow where the concept of being "at rest" 
is absolutely defined. Indeed, lattice models of multi-lane traffic flow have 
recently been studied [14]. The present model could be interpreted as a mul- 
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tilane flow in which vehicles cannot pass one another without an empty lane 
being present between them as they do so. However, the rules adopted for 
determining cluster velocities would themselves then have to be chosen non- 
Galilean- invariant (a cluster would travel only as fast as its slowest member). 
Apart from rotation effects [11], the presence or absence of Galilean invariance 
is (arguably) the major difference between colloidal and vehicular jamming. 

In summary, the counter scheme [10] is a deterministic way of implement- 
ing cluster flow which circumvents some of the difficulties associated with 
a stochastic implementation of cluster flow. In addition, it only allows single 
step movements of clusters and thus avoids potential problems associated with 
collisions of clusters and double occupancy of sites. 



5 Discussion 



In this work we have studied a simple stochastic model of shear-induced clus- 
tering on a 2d lattice. Our results strongly suggest that a discontinuous transi- 
tion occurs, with strong hysteresis, from an unjammed state to a jammed state 
as the shear rate is increased or temperature decreased. As the system jams, 
the distribution of clusters becomes power law. However, the exponent may 
have a non-trivial temperature dependence, appearing to vary continuously. 

Our conclusions are based on the study of finite systems and, while we have 
presented strong evidence in favour of the presence of a discontinuous phase 
transition in the model, care is required in identifying such a phase transition 
[7]. It would be useful therefore to study systematically finite-size effects in 
the system, as well as the effect of varying the aspect ratio. Also, we have 
not considered rotation of clusters but it should be possible to do so following 
ref. [11]. 

It is interesting to compare the present transition to several other models 
of reversible cluster-cluster aggregation in which the aggregation is due to 
diffusion of clusters, rather than an imposed shear. In the conserved-mass 
aggregation model of Krishnamurthy et al [15] a transition from a disordered 
(unjammed) phase to a phase containing an infinite cluster occurs as the rate 
of single particle dissociation increases. The infinite aggregate phase exhibits a 
power law cluster size distribution. However, in contrast to the present model, 
the transition in that work is reminiscent of Bose condensation and does not 
exhibit hysteresis. We also mention the work of Shih et al [16] and Haw et al 
[17] where diffusion limited cluster-cluster aggregation with finite bond energy 
(thus allowing fragmentation) is studied in the context of colloidal gels. 
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